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SUMMARY 

A simple matrix polynomial approach is introduced for approximating unsteady aerodynamics in the s-plane 
and ultimately, after combining matrix polynomial coefficients with matrices defining the structure, a matrix 
polynomial of the flutter equations of motion (EOM) is formed. A technique of recasting the matrix- 
polynomial form of the flutter EOM into a first order form is also presented that can be used to determine the 
eigenvalues near the origin and everywhere on the complex plane. An aeroservoelastic (ASE) EOM have 
been generalized to include the gust terms on the right-hand side. The reasons for developing the new matrix 
polynomial approach are also presented, which are the following: first, the “workhorse” methods such as the 
NASTRAN flutter analysis lack the capability to consistently find roots near the origin, along the real axis or 
accurately find roots farther away from the imaginary axis of the complex plane; and, second, the existing s- 
plane methods, such as the Roger’s s-plane approximation method as implemented in 1SAC, do not always 
give suitable fits of some tabular data of the unsteady aerodynamics. A method available in MATLAB is 
introduced that will accurately fit generalized aerodynamic force (GAF) coefficients in a tabular data form 
into the coefficients of a matrix polynomial form. The root-locus results from the NASTRAN pknl flutter 
analysis, the ISAC-Roger’s s-plane method and the present matrix polynomial method are presented and 
compared for accuracy and for the number and locations of roots. 

ACRONYMS 

ASE = Aeroservoelastic, 

eig = MATLAB function to solve generalized eigenvalue problems, 

EOM = Equations of motion, 

GAF = Generalized aerodynamic forces coefficient, 

GM, GD and GK = Generalized mass, damping and stiffness matrices, respectively, 

HSR = High Speed Research program, 

invfreqs = MATLAB function for fitting transfer functions from frequency response test data, 

1SAC = Interaction between flexible Structures, unsteady Aerodynamics and active Controls software 
utility, 
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LSE = Least squared error, 

MATLAB = Matrix computational utility, 

NASTRAN = Finite element structural modeling and flutter analysis software utility, 
TCA = Technical Concept Aircraft, and 

ZAERO = Unsteady aerodynamic and flutter analysis software utility. 


1.0 INTRODUCTION AND BACKGROUND 


1.1. Motivation and Study of Uncharacteristic Aeroelastic-Analysis Results 

No doubt, frustration helped to develop the approach presented in this paper. The initial concern that led to 
the subsequent frustrations was an unusual and uncharacteristic result obtained from a flutter analysis of a 
proprietary-design aircraft. In performing a NASTRAN 1 ® flutter analysis [1] and assessing the resulting root- 
locus plot, a structural root in the negative half of the complex plane took an unusual path from its zero 
dynamic pressure location near the imaginary axis to its high dynamic-pressure location very close to a rigid- 
body mode. When the two modes came in close proximity of each other and showed no apparent aeroelastic 
interactions normally appearing under those circumstances, a study was conducted to understand the 
uncharacteristic aeroelastic results. After resolving this first concern, later another concern materialized by 
using an s-plane analysis. This second concern revolved around aeroelastic modeling issues and flutter 
analysis methods, which made it necessary to create a third new method in an attempt to confirm the original 
s-plane method results and also to try to resolve the questionable results of the NASTRAN flutter analysis. 

An s-plane methodology was employed to better understand the root-locus results, since experience has shown 
that s-plane methods, such as the Roger’s s-plane approximation as implemented in IS AC [2], (or more 
concisely, the ISAC -Roger’s method), do a better job of capturing the behavior of the roots in the complex 
plane near and around the origin than does a frequency-based flutter analysis method. The ISAC version of 
Roger’s approximation that was employed here incoiporates enhancements specific for flight control 
applications and does not reflect the original form of the Roger’s fit algorithm. After a long struggle to obtain 
the s-plane fit coefficients, which were marginally adequate, the s-plane flutter equations of motion (EOM) 
were generated and the eigenvalues were computed to obtain the s-plane root-locus results. Although the 
flutter dynamic pressure, frequency and the path of the flutter root matched those of the NASTRAN analysis, 
the root-locus paths of the other modes did not follow the NASTRAN results. The uncharacteristic path of the 
NASTRAN structural root moving close to the rigid-body mode did not take place in the s-plane analysis, 
which helped to alleviate the initial concern. However, a result that looked very puzzling in the s-plane 
analysis was a flutter root involving a slightly higher structural mode going unstable at a slightly higher 


1 NASTRAN “ is a registered trademark of NASA and manufactured by MSC. Software Corp., Santa Ana, CA, USA. 

The use of trademarks or names of manufacturers in this report is for accurate reporting and does not 
constitute an official endorsement, either expressed or implied, of such products or manufacturers by the 
National Aeronautics and Space Administration. 
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dynamic pressure, not seen in the NASTRAN analysis. Now the concern led to aeroelastic modeling issues 
and issues with the flutter-analysis methods. 

In assessing the two results, more questions surfaced, such as, is the modeling sufficient, is the s-plane fit 
sufficiently accurate, and are there enough generalized aerodynamic force coefficients computed, especially at 
the lower reduced-frequency range? The NASTRAN and the s-plane flutter analyses were re-run with about 
30% more reduced frequencies in each. After assessing the new results, again the same respective root-locus 
characteristics were obtained by both methods. With the rather large differences in the results obtained by the 
NASTRAN and s-plane analysis methods, there were still lingering concerns that the s-plane fits lacked 
sufficient accuracy to give completely believable s-plane flutter results. This worry stems from the fact that 
these s-plane coefficients were computed by a legacy code, which is cumbersome at a minimum and has not 
been updated or recompiled as operating systems, libraries and memory allocation strategies changed with 
time. By using this code, significant doubt was cast on the correctness of the results. 

Since the IS AC -Roger’s method does not consistently give good fits, a subsequent effort was pursued trying 
various other forms of s-plane approximations to obtain better fits. Finally, a polynomial form was developed, 
implemented and found to give quite adequate fits of any GAF at nearly all ranges of reduced frequencies, if 
the order of the polynomial is sufficiently high. This polynomial approach of approximating the GAFs in the 
s-plane required the development of a whole new set of codes to generate the first-order equations of motion 
and to perform the root-locus flutter analysis. 

The root-locus results were nearly identical for case of (1) the polynomial method with the almost exact s- 
plane fits and (2) the ISAC-Roger’s method with the marginal fits. In addition, the dynamic pressures for the 
two flutter points for each case were nearly the same. With the agreement of the two s-plane-method results 
and NASTRAN giving the uncharacteristic root-locus result, the reasons for the second concern were resolved 
as being an aeroelastic modeling issue. Clearly, the various results of this aircraft configuration and the 
exercising the various flutter-analysis methods brought to light and demonstrated issues of modeling the 
aeroelastic systems and the issues involving the limitations or the benefits of the various methods in 
performing flutter analysis. 

Unfortunately no flutter results of the proprietary-design aircraft configuration described earlier can be 
presented here. Instead, the High Speed Research (HSR) program’s Technical Concept Aircraft (TCA) aircraft 
configuration was used to exercise the various methods presented in the paper. The TCA configuration 
presented in the Numerical Results section of the paper exhibits many of the same aeroelastic issues as the 
proprietary-design aircraft, but, more importantly, the TCA configuration is available for open discussion. 
Unlike the problems that surfaced with the proprietary configuration, the results from the NASTRAN and s- 
plane methods using the TCA gave nearly the same flutter points and root-locus paths. Thus, the effort of this 
paper became more of a validation task to show how well the different methods compared, as well as, to show 
some of the unique differences and characteristics of the methods. Hopefully, by this means, the effort will 
demonstrate that the new method is another, perhaps, a more precise and accurate, tool to perform flutter 
analysis in the s-plane. 

1.2. Benefits and Limitations of Frequency-Based Methods 

When performing flutter and gust-response analyses on aeroelastic systems with unsteady aerodynamics, the 
analysis capabilities in both ZAERO 2 ® [3] and NASTRAN use frequency-based methods. The frequency- 
based methods, such as the doublet-lattice or the ZONA aerodynamic codes, are simple to use and robust in 

2 ZAERO® registered trademark of ZONA® Technology, Inc., Scottsdale, AZ 85258-4578. 
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generating accurate representations of linear unsteady aerodynamics. Simple trapezoidal boxes, which these 
methods apply to form their aerodynamic grid patterns, can be conformed to a variety of very complex aircraft 
configurations with many non-planar surfaces and lifting surfaces oriented in many different directions. The 
convenience and versatility of the frequency-based methods become apparent when one can obtain flutter 
solutions of very acceptable accuracy by using simply generated unsteady aerodynamics, generalized masses 
and frequencies. 

There are some significant limitations of the frequency-based aeroelastic methods, such as their inability to 
find roots near the real axis or in the proximity of the origin of the complex plane, which is the case when 
investigating “body-freedom” flutter or static instabilities, such as divergence. This particular difficulty stems 
from the formulation of the basic eigenvalue problem itself for computing flutter where the frequency- 
dependent unsteady aerodynamic forces are available only for the positive half of the jco axis. Normally, 
when studying the roots located at or near the real axis to obtain representative root-locus results, it is 
important to model also the negative half of the j CO axis. The NASTRAN flutter analysis does not do this. 

1.3. Existing s-plane Methods and Their Limitations 

Casting frequency-dependent aerodynamics into the s-plane allows: the representation of the unsteady 
aerodynamics in the entire complex plane, the use of modem state -space methods to conveniently deal with 
aeroservoelastic systems; and the efficient and precise extraction of the eigenvalues on or near the origin, 
along the real axis and regions on and farther away from jco axis of the complex plane. Many of the s-plane 
methods such as the ones developed by Roger [4] and the minimum-state method developed by Karpel [5] 
generally do a good job of representing the frequency-dependent aerodynamics, but these methods are not 
straightforward to use, are labor intensive and require many iterations with fit parameters to obtain good fits. 
Many times, good or even acceptable fits are difficult to obtain or, in some cases, non-existent for complex 
aircraft configurations and mode shapes. The difficulty is compounded when attempting to fit the unsteady 
aerodynamics of the gust generalized forces. The need to develop a more efficient and accurate way of 
representing the unsteady aerodynamics in the s-plane and to harness the capabilities of modem and more 
efficient state-space methods inspired the methodology development presented here. 

1.4. Analytic Functions and Their Associated Properties 

Complex-variable theory [6] provides some insights of why the polynomial s-plane approximation produces 
the extraordinarily good fits of the unsteady aerodynamics. A very basic and crucial property that applies for 
GAF tabular data is whether the data can be represented as an analytic function. In complex-variable theory, 
the definition of an analytic function is that it must be continuous and its first derivative must also be 
continuous. Complex-variable theory further shows that, if a function is analytic, it can be represented by an 
infinite power series. Polynomials, in the context of the application presented here, represent a tmncated 
power series. This analytic -function property is very useful, since all valid tabular GAF data falls into this 
category. 

1.5. Goodness-of-Fit Criterion 

It is possible to define a criterion for goodness of fit, such as a not-to-exceed value of the least squared error 
(LSE) between the tabular values and the fitted values of the GAFs. In some s-plane fit procedures employed 
at Langley Research Center the LSE is computed, and it could be computed for the matrix polynomial 
method. However, in practice the LSE has been used as a criterion only rarely. For decades the standard 
practice has been for the engineers to visually examine the fits to see how well they match the tabular values, 
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regardless of the particular value of LSE. Engineering judgment has evolved into the primary (if not the sole) 
determining factor for the goodness of s-plane fit. 

2.0 DEVELOPMENT OF THE MATRIX POLYNOMIAL APPROACH 
2.1. GAFs in Polynomial Form 

In this paper, equation (1) is the new polynomial form proposed to represent the generalized aerodynamic 
forces (GAFs) in the s-plane. This form has been found to accurately “fit” the frequency-dependent GAFs, if 
the order of the polynomial is sufficiently high, 

GAF( k h ) ~ GAF(s ) (la) 

GAF (s ) = (A n s n + ... + A 2 s 2 + AjS + A 0 ) (lb) 



where, 


GAF (k h ) = tabular form of the matrix of the generalized aerodynamic forces, 

k h = tabular reduced frequencies, 
h = index for the tabular reduced frequency, 

GAF (s ) = polynomial approximation of the matrix of the generalized aerodynamic forces, 
Aj = matrix of the polynomial coefficients of the i th power of 5 , 

5 = reduced-frequency-scaled Laplace variable that is equal to T r s, 
s = Laplace variable scaled to the physical frequencies, 



V= vehicle velocity, and 
C r = reference length. 

2.2. Polynomial Fit Procedure 

The numerical procedure used to fit the coefficients for each element in the GAF matrix employs 
MATLAB 3 ®’s [7] invfreqs function that uses a fit procedure developed by Levi [8]. The invfreqs function 
converts easily and in a straightforward fashion the tabular magnitude and phase data into a transfer function 
of the form, 


’MATLAB® Registered trademark of The MathWorks, Inc, Natick, Massachusetts. 
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ORGANIZATION 


h( s) = 


As"+A n _ 1 s"- I +-- + A 0 
B d s d +B d y- 1 + --- + B 0 


(2a) 


However, for the application presented in this paper, only the numerator coefficients are obtained. All 
denominator coefficients are set identically to zero, except for the term B„ which is set to one. 

This simplification reduces the equation (2a) to the following form: 

h(s) = A n s n +A n _ lS nI + --- + A 0 (2b) 

In preparation for implementing the fit procedure for approximating GAFs, a change of variable from .s' to .S' is 
necessary, yielding the form, 


h(s) = A n r+A n _ 1 s n ~ I + --- + A 0 (2c) 

The invfreqs function employs a least squares procedure to solve directly for a minimizing set of the 
coefficients according to, 


in 


n 


min 

A 


I 

h=l 


GAF(k,)-'£ 


1=0 



(3) 


where, 


* = A- 

There are no fit constraints set in the way the coefficients are computed, except for the user to define the 
highest order of the polynomial used to perform the fits. 

2.3. Flutter Equations Formulation 

As a reference, the flutter EOM expressed in terms of the Laplace variable is presented first, 

[GMs 2 + GDs + GK + qGAF (T r s )] x = 0 (4a) 

where, 

q = l / 2 pV 2 , 

p = atmospheric density, 

GM = generalized mass matrix, 
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GD = generalized damping matrix, 

GK = generalized stiffness matrix, and 

x = generalized coordinate vector. 



If for example, a 4 th order form for the equation (2c) is substituted for the GAF(T,s) in equation (4a), the 
flutter EOM becomes, 

[qT 4 A 4 s 4 + qT 3 A 3 s 3 + (qT 2 A 2 + Gm)s 2 + (qT r A, + GD)s + (qA 0 + GK )]x = 0 (4b) 

In the equation (4b), the terms, [q T 2 A 2 + GM) , (qT r A, +GD) and (qA g +GK) are referred to as the 
augmented generalized mass, damping and stiffness matrices, respectively. 

Equation (4b) was recast into a first-order equation form to allow using a general-purpose eigenvalue 
computational procedure for the efficient solution of the eigenvalues. Additionally, a special canonical form 
of the first-order equations was chosen for this application because it has the property of being numerically 
well-conditioned when evaluating eigenvalues of large matrices, 


AEX = AX (5) 

where, X and X are a scalar eigenvalue and its corresponding eigenvector, respectively, and where, 

qT 4 r A 4 0 0 0 

E - 0 100 

0 0 10 

0 0 0 1 

-qT 3 A 3 I 0 

A = -{qT 2 A 2 + GM ) 0 I 

- {qT r A 1 + GD ) 0 0 

- (q~A a + GK ) 0 0 

To solve the generalized eigenvalue problem of equation (5) for the complex roots, MATLAB’s eig function 
was used. This function is specifically set up to solve eigenvalues of state-space systems. 

As seen by examining equations (5), (6) and (7), the scaling parameter, T r , plays the important role of scaling 
the EOM to the physical frequencies. The scaling parameter raised to the appropriate power multiplies each 
of the polynomial coefficients matrices of the s-plane fit. Because of V residing in the denominator of T„ T r 
raised to progressively higher integer powers, T" , results in progressively smaller values of the quantity, T" , 


0 

0 

I 

0 


( 6 ) 


(V) 
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and thereby progressively diminishes the influence of the higher order polynomial coefficient matrices. This 
effect is amplified as the airspeed is increased. This is a fortuitous result because it makes the fit procedure 
especially attractive in that, potentially, a smaller number of the lower order aerodynamic terms may be 
sufficient to accurately represent the aeroelastic system. At this point of the methodology development, the 
highest order needed for sufficiently accurate representation of aeroelastic system has not completely been 
determined. 

2.4. Aeroservoelastic Equations Formulation 

A 4 th order form of the aeroservoelastic (ASE) equations is presented in equation (8). The left-hand sides of 
equation (4a) and (8) are identical. The right-hand side of equation (8) contains the control and gust terms. 

[qT 4 A 4 s 4 + qT 3 A 3 s 3 + (qT 2 A 2 + Gm)s 2 + (qT r A, + GD)s + (qA n + GK)]x = 

- R T r A c/ + (t T rA/ + {H T r A c2 + GM c V + (<l T A 1 V + (<lAo )J X c (8) 

- R T r A 8 4 s4 + 9 r r A g3 s 3 + (qT 2 A g2 )s 2 + (qT r A gl )s + (qA g0 )]w g 



where, 


A ci = matrix of polynomial coefficients of the i th power of s for the control term, 

A g i = matrix of polynomial coefficients of the i th power of s for the gust term and 
GM C = matrix of the control-surface generalized inertias and the mass coupling terms. 


The A c and A g matrices have the same number of rows as the coefficients on the left side of the equation, but 
the number of columns of A c depends on the number of control inputs, x c . The gust input is generally the 
vertical component of the atmospheric turbulence, w g . 

For numerical simulations and time-domain analyses purposes, the following first-order form of the 
aeroservoelastic (ASE) equations is presented, which is more useful than the polynomial form of equation (8), 


( 9 ) 


The E and A matrices in equation (9) are those given in equations (6) and (7). The additional matrices for the 
control and gust inputs are as follows: 


Ex = Ax + B c 

u c 

+ B 

V 


A. 

g 

_*v 


B 


qT 3 A c3 qT 4 A c4 

(qT 2 A c2 + GM c ) 0 
qT r A cI 0 

qAo o 


( 10 ) 
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0 

This formulation of the ASE equations only requires x c and dx/dt and w g and dw/dt inputs, which can be a 
significant advantage compared to other formulations. When the 1SAC -Roger’s s-plane fits are employed, the 
resulting formulation of the ASE equations [9] also requires d 2 ( )/dt 2 inputs for both the control and gust, 
which are acceleration variables that may not always be available in simulations. Thus, the lack of these 
acceleration variables can make the modeling of ASE input dynamics incomplete. These particular terms are 
especially important when applied to the control surfaces, since they pass on the control-surface inertial forces 
to the ASE model. 

3.0 NUMERICAL RESULTS 


and B 


« T r A g s 

9T r 2 A g2 

9 T r A gl 

. Ho 



The analytical model used in this study 
is based on the Technical Concept 
Aircraft (TCA) (Figure 1), which is one 
of the configurations investigated during 
the NASA High Speed Research 
program. In 1997, Boeing-Long Beach 
created an MSC NASTRAN finite 
element-model under a contract to 
NASA Langley Research Center. It 
represents a free-free, strength-sized 
aircraft configuration developed to study 
several aeroelastic issues including 
body-freedom flutter (BFF). With a 
great deal of effort by Boeing, much of 
the BFF behavior was eliminated to a 
level considered acceptable within its 
flight envelope. However, some residual 
tendencies toward BFF remain in the analytical model. For this reason, the TCA configuration was deemed to 
be a good choice for an evaluation of the matrix polynomial approach for approximating unsteady 
aerodynamics. 

3.1. NASTRAN Flutter Analysis Results 

The plots in figure 2 show the root locus results from an altitude-variation NASTRAN pknl flutter analysis 
using 25 structural modes. In the altitude variation, the altitude parameters at Mach 0.80 were varied from a 
very high altitude where the air density is nearly zero to a density equivalent to 40,000 feet below sea level 
that generated the maximum dynamic pressure of 3,300 psf. The altitude parameters for the flutter analysis 
were selected to vary the dynamic pressure from a range where the aerodynamics has nearly zero effect to 
where it has a very large effect. The scales for the real and the imaginary axes in the figure 2(a) were chosen 
such that the root loci of the entire complex plane can be directly compared in part (a) of figures 2, 4, 6 and 8 
resulting from the NASTRAN, 1SAC implemented ISAC-Roger’s method and the matrix polynomial method. 
The scales in figure 2(b) were chosen such that the root loci in part (b) of figures 2, 4, 6 and 8 can be directly 
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compared near the complex-plane origin resulting from each method. In figure 2(b), the first fuselage bending 
mode, whose root-locus starts at 7.0 on the imaginary axis, couples with the short-period mode as indicated by 
its trace curving down and then back toward the imaginary axis. Flutter is seen to occur at a dynamic pressure 
of 640 psf as the second structural mode crosses into the right-half plane. The iterative method uses, in its 
flutter computations, only the frequency dependent aerodynamics defined in the imaginary axis for its 
estimation of the aerodynamics throughout the complex-plane resulting in imprecise aerodynamics and 
therefore inexact root locations, especially those located farther away from the jco axis. For this flutter 
analysis, NASTRAN tracks the short period, which is a rigid-body mode, very well, generally not expected of 
a method using a frequency-based flutter root finding technique. Not withstanding, it should be 
reemphasized, NASTRAN does a very good and consistent job of finding structural roots on or in close 
proximity to the imaginary axis, which is the area most important in finding flutter roots and the flutter 
conditions. 

3.2. ISAC implemented Roger’s s-Plane Method Results 

The fits using ISAC-Roger’s s-plane method with four lag terms is shown in figure 3 for GAF(4,3) and 
GAF(4,5) from displacement mode 4 and pressure modes 3 and 5, respectively. There appears to be a slight 
discrepancy between the tabular values and the fitted results, but such differences are expected and are 
normally deemed well within the acceptable range for the ISAC-Roger’s method. In general, using ISAC- 
Roger’s s-plane method, it may take many tries with many adjustments of the fit procedure parameters to 
obtain acceptable fits. 

As with the NASTRAN flutter analysis, an altitude-variation flutter analysis was performed and the root-locus 
results at Mach 0.8 are presented in the plots of figure 4. Figure 4(a) shows the root-locus results using ISAC- 
Roger’s method, including the complex-conjugate roots in the bottom half of the complex plane. The 
difference between the ISAC-Roger’s method root-locus results and the NASTRAN root-locus results may be 
caused by the imprecise method that the NASTRAN pknl flutter analysis uses to estimate unsteady 
aerodynamic forces for computing roots farther away from the imaginary axes. Figure 4(b) contains the first 
eight roots and can be compared to figure 2(b) near the complex plane origin. As indicated in the figure, the 
flutter-point dynamic pressure is shown at 620 psf, about 3% lower than the NASTRAN’s result. In 
comparing figure 4(b) to figure 2(b), it can seen that, in NASTRAN, the roots appear to travel much farther to 
the left in the complex plane than those of the ISAC-Roger’s s-plane method. In this plot, as in the 
NASTRAN results, the coalescing between the first fuselage bending mode, labeled as BFF, and the short- 
period modes is readily apparent as the root-locus is also taking a downward curved track toward the 
imaginary axis. Of course, the path from the origin of the short-period mode differs slightly by not moving 
quite as far up the imaginary axis or as far out the negative real axis. 

3.3. Matrix Polynomial Approach Results 

This section of the paper presents results obtained using the matrix polynomial approach. First a 20 th -order 
polynomial fit is employed. Then a 3 rd -order polynomial truncation of the 20 th -order fit is employed. 

Figure 5 contains the 20 th -order polynomial fits of GAF(4,3) and GAF(4,5). In comparing the GAF(4,3) and 
GAF(4,5) fits presented in figure 3 to those of figure 5, it is clear that the polynomial fits are much closer to 
the tabular points then are the ISAC-Roger’s fits. A closer examination of the real and imaginary GAF(4,5) 
plots of figure 5 reveals that the fit looks nearly exact. However, both the real and imaginary parts of the 
GAF(4,3) plot shows some slight deviation of the polynomial function between the tabular points of the 
highest two reduced frequencies. Surveying the plots of nearly all the GAFs, similar, very small deviation 
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results were obtained at the highest reduced frequencies for some GAFs. 

Although not shown, the resulting 20 th -order flutter EOM similar to equation (4b) was put into a canonical 
first-order form similar to equations (6) and (7) so that the flutter solutions could be obtained by solving the 
generalized eigenvalue problem, equation (5). The resulting order of the generalized eigenvalue problem, 
and, therefore, the number of roots in the solution, was 440. This compares to 25 for the NASTRAN solution 
and 125 for the IS AC -Roger’s method solution. 

Figure 6 contains the corresponding flutter root locus. Figure 6(a) may be compared to figure 2(a) and 4(a). 
The most obvious feature of figure 6(a) is the large number of roots, about half of them unstable, associated 
with the 20 th -order s-plane fits. The vast majority of these roots occur in complex-conjugate pairs; the rest are 
located on the real axes. Although, the method represents the unsteady aerodynamics in the complex plane 
very well, when the flutter equations are solved, besides finding the roots important to determine flutter points 
and conditions, the eigensolver finds all the roots including a very large numbers of aerodynamic roots. In 
most cases, these roots are situated generally well outside the area of interest for flutter and can easily be 
distinguished and ignored. 

Figure 6(b) may be compared to figure 2(b) and 4(b). In overall appearance the root locus in figure 6(b) is 
very similar to the other two root loci. The interaction of the flexible mode and the short-period mode is the 
same and the flutter dynamic pressure is the same as that predicted using ISAC-Roger’s fits. 

The matrix polynomial approach imposes no constraints regarding the stability of the aerodynamic roots. 
However, methods such as ISAC-Roger’s generally produce stable aerodynamic roots by the constraints 
imposed during the fit procedure. In dealing with the large number of aerodynamic roots produced by the 
matrix polynomial approach, there are a few computational options available to reduce size of the equations of 
motion to abrogate the effects of the aerodynamic roots. These techniques involve using model residualizing 
such as Guyan reduction methods; modal reduction methods, which eliminates unwanted eigenvalues; using 
iterative root-finding techniques such as the pknl method of NASTRAN; or simply truncating the higher order 
polynomial coefficients. 

To demonstrate the polynomial coefficient truncation technique, a third-order polynomial formulation of the 
flutter equations of motion is given as an example. Previous experience of the author has shown that starting 
with a third-order polynomial fit to produce a third-order flutter EOM generates poor flutter results. To obtain 
relatively good flutter results, it is necessary to fit a full order polynomial, for this example a 20 th -order, then 
use only the lowest order terms for the third-order flutter EOM. Figure 5 shows the results of the full order 
polynomial form of the GAF using all the coefficients, and figure 7 shows the results using the same set of 
coefficients, but applying only the third-order terms and lower. At the lowest reduced frequency range, the 
20 th -order fit truncated to third-order polynomial generally follow the tabular GAFs but, at the higher range of 
reduced frequencies, the truncated polynomial results substantially deviate from the tabular results. 

In comparing figure 8(a) and figure 6(a), the number of aerodynamic roots produced by the third-order flutter 
EOM is considerably less than those produced by the 20 th -order EOM. Also the only unstable aerodynamic 
roots present in figure 8(a) are a cluster along the positive real axes. This outcome provides a more 
compelling incentive to reduce the polynomial order for both the flutter and ASE EOM. Now comparing 
figures 8(b) and 6(b), again the appearance the root locus are similar, but the flutter dynamic pressure is 
substantially lower, from 620 psf for the 20 th -order case to 560 psf for the third-order case. Near the origin, 
the rigid-body mode root locus of the third-order case matches the 20 th -order case. 
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Figure 9 shows a close-up of these rigid-body modes for the NASTRAN, the ISAC -Roger’s method, the 20 th - 
order and the third-order cases. In the inertial coordinate system, the rigid-body modes involve the pitch and 
plunge dynamics. The pitch mode is characterized by the short-period root-locus paths shown in figures 2(b), 
4(b), 6(b) and 8(b). A flutter analysis of a free-free aircraft configuration should normally depict the plunge 
mode as two roots near the origin on the real axes. This feature was captured neither by NASTRAN nor by 
the ISAC-Roger’s method. Using the pknl method of NASTRAN, the plunge root-locus never appears 
because of the un-modeled aerodynamics at zero reduced frequency and the iterative method used to find 
roots. For the ISAC-Roger’s method, the plunge roots stay at the origin because the plunge aerodynamic 
forces are constrained to be zero when fitting the plunge unsteady aerodynamics, even though these 
aerodynamic forces are not computed to be exactly zero. When fitting with polynomial functions, there are no 
constraints set for the fit of any of the polynomial GAFs including the polynomial GAF of the plunge 
aerodynamics. With the plunge aerodynamics included in the flutter EOM, a subtle but important difference 
appears between the root locus of the polynomial method and the other two methods. At very low dynamic 
pressure, the neutrally stable plunge rigid-body mode at the origin in figure 6(b) and figure 8(b) separate into 
two roots: one stable and one unstable. With increasing dynamic pressure, these two plunge roots are taking 
trajectories on the real axes in opposite directions away from the origin. The motion of the unstable plunge 
root may be visualized as a free-free aircraft deviating away from a reference plane and the increasing values 
of root-locus points farther to the right of the origin correspond to the increased exponential rates at which the 
aircraft departs from the reference plane. 

4.0 CONCLUSIONS 

A method using matrix polynomials has been developed to more accurately fit the unsteady aerodynamics in 
the flutter equations of motion to find eigenvalues at the origin and everywhere else on the complex plane. 
Through a numerical fit procedure, the method is found to almost exactly represent the s-plane using tabular 
generalized aerodynamic force (GAF) coefficient data generated by linear unsteady aerodynamic codes. From 
complex-variable theory, the polynomial equation is mathematically a general form and well suited to fit any 
tabular aerodynamics data. When compared to results of other s-plane fit methods, the matrix polynomial fit 
procedure is straightforward and easy to use. Fit results from the ISAC-Roger’s method, from a high-order 
polynomial fit and from a low-order truncation were presented to demonstrate the accuracy of the fit 
procedures. Flutter root-locus results have been presented. Root-locus results of rigid body modes were 
presented to show the increased fidelity attained by using the matrix polynomial approach to model the rigid- 
body dynamics. 
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Figure 2. Root locus plot using NASTRAN pknl Flutter Analysis at M = 0.8. 
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Figure 3. Typical fit of TCA unsteady aerodynamics at Mach = 0.8 using 
ISAC-Roger’s s-plane approximation. 
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Figure 4. Root locus plot of flutter analysis with s-plane using ISAC-Roger’s approximation 
aerodynamics at Mach = 0.8. 
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Figure 5. 20 th -order polynomial fits and tabular values of the GAF(4,3) 
and GAF(4,5) at Mach = 0.8 of a TCA configuration. 
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Figure 6. Root locus plot of flutter analysis with s-plane method using 20 th -order matrix polynomial 
approximation of the unsteady aerodynamics at M = 0.8. 
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Figure 7. 20 th -order fit truncated to third-order polynomial and the tabular 
values of the GAF(4,3) and GAF(4,5) at Mach = 0.8 of a TCA configuration. 
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Figure 8. Root locus plot of flutter analysis with s-plane method using third-order 
matrix polynomial approximation of the unsteady aerodynamics at M = 0.8. 
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Figure 9. Root locus plots of the rigid-body modes from various flutter analyses at M = 0.8. 
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